% fig3a.m ——— 复刻论文 Fig.3(a)，可自定义截断与绘图区间 ———
clear; clc; close all;

%% —— 用户自定义参数 —— 
timeDiscard  = 200;   % 丢弃前 200 单位时间
timePlotBeg  = 200;   % 绘图起始时间
timePlotEnd  = 500;  % 绘图结束时间

%% —— 模型参数 —— 
alpha = 3;
beta  = 1;
xi    = 0.8;
mu    = 1;
A     = 3.4;
omega = 1;

%% —— 仿真设置 —— 
dt  = 0.01;             % 时间步长
T   = timePlotEnd;      % 总仿真时间至少要 >= timePlotEnd
N   = floor(T/dt);
t   = (0:N-1) * dt;     % 时间向量

% 计算各个时间点对应的下标
idxDiscard = floor(timeDiscard  / dt) + 1;
idxBeg     = floor(timePlotBeg  / dt) + 1;
idxEnd     = min(floor(timePlotEnd / dt) + 1, N);

%% —— 初始化 —— 
S = zeros(3, N);
S(:,1) = [0.2; 0.1; 0.2];
H = zeros(1, N);

%% —— RK4 仿真 —— 
for k = 1:(N-1)
    mu_s = A * cos(omega * t(k));
    %—— 主系统 RK4 —— 
    k1 = f(S(:,k),        alpha,beta,xi,mu,mu_s);
    k2 = f(S(:,k)+0.5*dt*k1, alpha,beta,xi,mu,mu_s);
    k3 = f(S(:,k)+0.5*dt*k2, alpha,beta,xi,mu,mu_s);
    k4 = f(S(:,k)+    dt*k3, alpha,beta,xi,mu,mu_s);
    S(:,k+1) = S(:,k) + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
    %—— 瞬时能量 —— 
    xk = S(1,k);  yk = S(2,k);
    H(k) = xk^2/(2*alpha) + yk^2/2;
end
% 末点能量
xk = S(1,end); yk = S(2,end);
H(end) = xk^2/(2*alpha) + yk^2/2;

%% —— 截断并提取绘图数据 —— 
t_plot = t(idxBeg:idxEnd);
Sx = S(1, idxBeg:idxEnd);
Sy = S(2, idxBeg:idxEnd);
H_plot = H(idxBeg:idxEnd);

% 如果要去掉丢弃段对平均能量的影响，可这样计算：
Hmean = mean(H(idxDiscard+1:idxEnd));

%% —— 绘图 —— 
figure;

% (a1) 相位图
subplot(1,3,1);
plot(Sx, Sy, 'b-', 'LineWidth',1);
xlabel('x','FontSize',12); ylabel('y','FontSize',12);
title('(a1) 相位图','FontSize',14);
grid on;

% (a2) x-τ 时序
subplot(1,3,2);
plot(t_plot, Sx, 'b-', 'LineWidth',1);
xlabel('\tau','FontSize',12); ylabel('x','FontSize',12);
title('(a2) x-\tau 时序','FontSize',14);
grid on;

% (a3) 能量 H(τ)
subplot(1,3,3);
plot(t_plot, H_plot, 'b-', 'LineWidth',1);
xlabel('\tau','FontSize',12); ylabel('H','FontSize',12);
title('(a3) 能量 H(\tau)','FontSize',14);
grid on;
% 标注平均能量
text( t_plot(end)*0.3, max(H_plot)*0.8, ...
      sprintf('\\langleH\\rangle=%.3f', Hmean), ...
      'FontSize',12, 'Color','k');
